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The properties of satellite galaxies in simulations of galaxy 
formation 
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ABSTRACT 

We investigate the properties of satellite galaxies in cosmological iV-body/SPH simula- 
tions of galaxy formation in Milky Way-sized haloes. Because of their shallow potential 
wells, satellite galaxies are very sensitive to heating processes which affect their gas 
content. Their properties can therefore be used to constrain the nature of feedback 
processes that regulate galaxy formation. In our simulations, we assume that all the 
energy produced by supernovae is used as kinetic energy to drive galactic winds. Sev- 
eral of our simulations produce bright, disc-dominated galaxies. We find that wind 
models in which the wind speed, v w , is proportional to the local velocity dispersion of 
the dark matter, a (and thus the wind mass-loading, rj w cx c~ 2 ), make star formation 
in satellites sporadic, reproduce the observed satellite luminosity function reasonably 
well (down to M v — —7) and match the luminosity-metallicity relation observed in 
the Local Group satellites. By contrast, models that assume a constant wind speed 
overproduce faint satellites and predict an incorrect luminosity-metallicity relation. 
Our simulations therefore suggest that the feedback processes that operate on the 
scale of satellite galaxies should generate galactic outflows whose mass-loading varies 
inversely with the depth of the potential. 

Key words: methods: numerical - galaxies: evolution - galaxies: formation - cos- 
mology: theory. 



1 INTRODUCTION 

The 'A-cold dark matter' (ACDM) model is now broadly ac- 
cepted as the standard paradigm in cosmology. This model 
completely specifies the initial conditions for the formation 
of cosmic structure, as well as the values of the cosmological 
parameters. Thus, in principle, it is possible to calculate the 
predicted evolution of objects of any kind, provided that the 
relevant astrophysical processes are understood. The com- 
plexity of these processes makes computer simulations the 
ideal tool to investigate the formation of structure in the 
Universe. 

There is a long history of attempts to simulate the 
formation of spiral galaxies from ACDM initial con- 
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on an 'angular momentum problem' was identified. In the 
ACDM model, galaxies are generally assembled by mergers. 
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As fragments merge, their angular momentum, which is 
primarily invested in their orbits, is transferred to the outer 
halo. Thus, if these fragments already carry within them the 
majority of the baryons which will make up the final galaxy, 
then these baryons will end up in a slowly rotating central 
spheroid rather than in a disc (Fre nk et al.l [l985). The 
likely solution to the problem was also identified early on: 
the baryons must be kept from cooling into the fragments, 
perhaps by some form of feedback, until the merging 
activity has subsided after which they can cool smoothly to 
form a disc (jNavarro fc Bendll99ll ; INavarro fc Whitd[l99l : 
INavarro. Frenk. fc Whitdll995h . 

The nature of the angular momentum problem and 
the key role of feedb ack were explicitly demonstrated by 
lOkamoto et~afl (|2005l ) who simulated the formation of three 
galaxies - all in the same dark matter halo - assuming 
three different models of feedback powered by the energy 
liberated in supernovae (SNe). The resulting three galaxies 
spanned the whole range of morphological types, from a pure 
sph eroid to a disc gala xy with bulge-to-total mass ratio of 
0.5. IZavala et al.l (|2008T ) explored these galaxies further and 
showed how, when the feedback is not strong enough to pre- 
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vent gas cooling early on, a spheroid forms, whereas strong 
feedback can result in the formation of a disc by infall of gas 
which approximately conserves its angular momentum. 

The physical processes that establish the feedback 
loops required for disc galaxies to form are poorly un- 
derstood, but almost certainly involve SNe energy and 
perhaps al so energy extracted from a central black 



hole (e.g. |Pi Matteo et all 120051; lOkamoto et all 120051: 



ISiiacki et all 120071 : iBooth fc Schavd 120091 ). These processes 
occur on scales which are many orders of magnitude be- 
low those that can be followed directly in galaxy simu- 
lations from ACDM initial conditions. They must there- 
fore be treated as 'sub-grid' physics. A number of au- 
thors have attempted this kind of modelling using var- 
ious recipes leading to a new generation of simulations 
which are producing increasing l y realistic disc galax- 



ies dSpringel fc Hernquistl l2003al: iRobertson et al.l 20041: 



Okamoto et all 120051: IStinson et alj|2006l; iGovernato et al l 



20071 ; IScannapieco et alll2008l : ICeverino fc Klvpinll2009h . 

Although a full physical understanding of how galactic 
feedback processes work seems a long way away, it is worth 
trying to identify or constrain the 'macroscopic' properties 
to which those processes must give rise for a galaxy simula- 
tion to be successful. Since it seems clear that much of the 
activity must take place early on, a promising approach is to 
investigate how different feedback models influence the prop- 
erties of surviving fragments which appear today as satellites 
of the main galaxy. This is the approach we take in this pa- 
per: we focus on the luminosity function and metallicity of 
satellites and examine how these properties are affected by 
different models of feedback. 

The properties of satellites are interesting in their own 
right. For example, their luminosity function has often been 
claimed to pose a major problem for the CDM model be- 
cause the number of satellites seen around the Milky Way 
is several orders of magnitude smaller than the number of 
substructures that s urvive the collapse of galactic haloes i n 
A-body simulations (|Klypin et al.ll 19991 ; iMoore et al.lll999h . 
However, there are various feedback processes that could 
explain why only a tiny fraction of subhalo es manage to 
make stars. SN heating is clearly one of them (| Benson et all 
2002, 2003); another is the photoionization of gas at early 
times, which suppresses gas cooling and raises the thermal 
pressure of the gas, preventing it from ac c reting onto small 



dark matter haloes iBa bul & Rccs 1992; Efstathiou 



Thoul fc Weinberg! 1 19961; iQuinn. Katz. fc E fstathiou 



1992; 



1996; 



Gnedinll2000l ; lOkamoto. Gao. fc TheundbOOSar T 

These processes can be readily modelled using semi- 
analytic techniques and such studies have shown that in- 
deed they greatly reduce the numbe r of visible satellites 
llBullock. Kravtsov. fc Weinberg! I2OO0I; iBenson et all |2002| ; 
ISomervillell2002l ; lLi et all 120091 ; iMaccio' et al1l2009h . Some 
doubt has recently been cast on the validity of these impor- 
tant conclusion because b oth the early st udie s and the more 
recent , but similar ones bv lLi et all (|2009l ) and lMaccio' et all 
(2009), were based on Gnedin's (2000) filtering mass de- 
scription of reionization which appears to overe stimate sig- 
nificantly the suppression o f galaxy formation faoeft et all 
120061 ; lOkamoto et al.ll2008al ). lMaccio' et all (|2009l ). however, 
have sho wn that even with the weaker suppression advo- 
cated by I Okamoto et all l|2008al ). acceptable satellite lumi- 
nosity functions can still be obtained. 



Hydrodynamic simulations are well suited for studying 
complex processes such as reionization and SN feedback ab 
initio, with a minimum of simplifying assumptions. An ear- 
lier generation of galaxy simulations showed that at least the 
brigh test satellite galaxies are produced in acceptab le num- 
bers (|Governato et al.l 120071 : iLibeskind et al.l 120071 '). These 
simulations, however, lacked the numerical resolution re- 
quired to study the bulk of the observed Milky Way satel- 
lites. 

The simulations presented in this paper include a num- 
ber of baryonic processes known to be relevant to galaxy for- 
mation. In particular, we model SN feedback as winds, a phe- 
nomenological approach which was originally introduced by 
ISpringel fc Hernquistl l|2003al ). IOp"penheimer fc Pavel (|2006h 
modified this treatment allowing the wind s peed and mass- 
loading to vary with galaxy properties, while lOkamoto et all 
l|2008bl ) introduced a model in which the energy in the winds 
reflects the timed-release of SN energy. Un like in previous 
simulations, Dalla Vecchia fc Schavel f20O8) considered the 
case in which wind particles are allowed to interact hydro- 
dynamically with interstellar medium (ISM) particles and 
found that due to hydrodynamical drag in the ISM, both 
the effective wind speed and the effective mass-loading (i.e. 
the wind speed and mass-loading when the wind leaves 
the star-forming region) change significantly from the input 
wind speed and mass-loading. This type of win d models are 
also used in recent cosmo logical simulations (e.g. lCrain et all 
120091 ; ISchave et al1l20ld ) . 

Winds are n ow observed in many star-forming galax- 
ies. Earlier data (|Martinl 1 19991 ; iHeckman et all |2000| ) sug- 
gested that their properties were indepen dent of galac- 
tic properties. More recent ob servations (|Martinl 120051 ; 
iRupke. Veilleux. fc Sanders! 12005). however, have reversed 
this view: galaxies with higher circular velocity and higher 
star formation rates appear to drive faster winds. In this pa- 
per, we consider two types of wind models, both of which as- 
sume that all the energy released from SNe is converted into 
the kinetic energy of the wind s. In one case, the wind speed i s 
assumed to be constant, as in ISpringel fc Hernquistl (|2003aD 
and as suggested by the earlier data ( Martin 19991 ). In the 
other case, the wind speed is assumed to be proportional to 
the lo cal velocity dispersion, as sugges ted by the more recent 
data (jMartin 2005; Rupke et al.ll2005l ). Similar, but different 
wind models have been used s uccessfully to investigate the 
eject ion of intergalactic metals (|Oppenheimer fc Davdl200e1 . 
2008) and the galaxy mass-metallicity relation of massive 
galaxies (M* > 1O 9 M ; iFinlator fc Dave1l2008l ). 

This paper is organised as follows. In Section 2, we de- 
scribe our simulations, providing brief descriptions of our 
modelling of baryonic processes (Section 2.1), the energetic 
and chemical feedback due to stellar evolution (Section 2.2) 
and the implementation of winds (Section 2.3). In Section 3 
we present our results. We first focus on the properties of the 
central galaxies (Section 3.1). We then address how the lumi- 
nosity functions and metallicities of the satellite galaxies are 
affected by different feedback prescriptions (Section 3.2). In 
Section 3.3 we perform convergence tests by using a higher 
resolution simulation. Finally, in Section 4, we discuss the 
results and summarise our main conclusions. 
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2 THE SIMULATIONS 

We selected the three most massive o bjects from the sampl e 
of six Aquarius haloes described in ISpringel et alJ l)2008t ) . 
The Aquarius haloes, extracted from a large cosmological 
simulation, were chosen to be relatively isolated and to have 
masses suitable to host a Milky Way-like galaxy. With one 
exception, which is not in our sample, they all have a rel- 
atively quiet merger history since 8 = 1. The three most 
massive haloes are: Aq-A, Aq-C, and Aq-D. Halo Aq-A 
has the quietest merger history of the whole sample and 
has been studied in the most detail to date. All our sim- 
ulations follow galaxy formation within a periodic cube of 
side 100 /i" 1 Mpc in a ACDM cosmology with parameters: 
fi m = 0.25, = 0.75, n b = 0.045, cr 8 = 0.9, n s = 1, and 
H = 100 h km s" 1 Mpc" 1 = 73 km s" 1 Mpc" 1 . These val- 
ues are the same as those a dopted for the Millennium Simu- 
lation l|Springel et al.ll2005h and the Aquarius project, and, 
within the uncertainties, are consistent with the current set 
of constrain ts from the WMAP 1 -year (|Spergel et alj|2003h 
and 5-year |Komatsu et alj|2009h data. 

Hydrodynamic simulations of galaxy formation in all six 
Aquar ius haloes have been presented by IScannapieco et all 
(2009). These have similar resolution to our simulations but 
use a different code and focus on different issues. They em- 
ployed a slightly lower value of the mean baryon density, 
fib = 0.04, than us. In that work, Aq-C, produced the most 
disc-dominated galaxy and had the smallest mass fraction 
in subhaloes within rso (the radius inside which the density 
is fifty times the critical density). Aq-D has a similar halo 
mass to Aq-A and Aq-CQ. The concentration parameter of 
this halo is the smallest of the three we have chosen Q. 

Since hydrodynamic simulations are much more com- 
putationally expensive than iV-body simulations with the 
same number of dark matter particles, we are only able to 
simulate ha loes at the lowest and next-to-lowest resolution 
employed in lSpringel et all (|200Sl ) (corresponding to levels 5 
and 4 in that paper). To test for numerical convergence, we 
carried out a high resolution simulation of Aq-D, to which 
we refer as Aq-D-HR (which corresponds to level 4). In Ta- 
ble [T] we summarise the halo properties, the particle masses 
and the gravitational softenings used for the simulations. 



2.1 Baryonic processes 

The simulation code was developed from an early version 
of GADGE T-3. This is similar to the publicly available 
GADGET-2 jSpringej^Oot ). but it has a more efficient load- 
balancing algorithm. We have modified GADGET-3 to follow 
a number of additional physical processes as follows. 

Both photo-heating by a spatially uniform, 
time-evolving ultraviolet background and radia- 
ti ve cooling depend on gas met allicity, as described 
in lWiersma. Schave. fc Smith! (|2009al ). by assuming the gas 
to be optically thin and in ionization equilibrium. We use the 



1 Other Aquarius haloes have mass a factor of ~ 2 smaller than 
these haloes 

2 The values of the concentration parameter, c, for haloes Aq- 
A, Aq-C, and Aq-D, obtained by fitting an NFW density pro- 
file (n avarro. Frenk, fc Whitelll997l) , are ~ 16, 15, and 9, respec- 
tivelv llSprineel et al.ll200g|) . 



ultraviolet background given by lHaardt fc Madau (12001 ) 
which is switched on at z = 9 (see also Okamoto et all 
(200 83)). The cooli ng/heating routine described by 
W iersma et al.1 (|2009al ) computes the cooling/heating rates 
individually for eleven elements whose contribution to the 
radiative cooling is significant (H, He, C, N, O, Ne, Mg, 
Si, S, Ca, and Fe). In the simulations, we track only nine 
ele ments and we take S a nd Ca to be proportional to Si 
fsee lWiersma et al.ll2009bh 

As in lOkamoto et~al] |2005l ). we use the smoothed 
metallicity, rather than the particle metallicity (which is 
normally used), to compute the photo- heating and radiative 
cooling rates and to give the initial amount of each individ- 
ual che mical element to newly-born stars. |w iersma et al.1 
(2009b) found that the particle metallicity underestimates 
the metallicity of relatively unenriched gas compared to the 
smoothed value and thus significantly underpredicts the to- 
tal stellar mass formed in their cosmological simulations. Us- 
ing the particle metallicity is particularly problematic when 
a small fraction of metal enriched particles are ejected as 
winds, as in our simulations. 

Sta r formation is modelled as in lOkamoto et al.1 
(2008b). That is, we treat an SPH particle whose den- 
sity exceeds a threshold density for star formation (tih > 
n,H,th = 0.1 cc" 1 ) as a hybrid particle that contains two 
distinct phases: hot ambien t gas and cold clouds, as in 
ISpringel fc Hern quist (2003a). In this model, the clouds have 
a fixed mass spectrum and the star formation rate in a gi- 
ant molecular cloud is a function of t he ambient pressure 
accord i ng to the two-pha se model of ISamland fc Gerhard! 
(20031). iBooth et al.1 (|2007i ) introduced an alternative to the 
hybrid particle approach in which the cold clouds are repre- 
sented by sticky particles (and hot, warm, and diffuse phases 
are represented by SPH particles). While this model shares 
many physical processes with our model, such as cloud for- 
mation by thermal instability and cooling and evaporation 
by thermal conduction, it allows the cloud mass spectrum 
to evolve with time through coagulation. 

Ste llar evolution is modelled as in lOkamoto et al.1 
(2008b]), with two exceptions: (i) we emplo y the more rea l- 
istic Chabrier initial mass f unction (IMF, Ichabrierl 120031 ) , 
instead of the Salpeter IMF (|Salpeter| [l955): and (ii) we use 
metallicity-dependent stellar lifetimes and chem ical yields 
j|Portinari. Chiosi. fc Bres san 19981; lMarigoll200ll ). The code 
does not assume the instantaneous recycling approxima- 
tion. The production of metals by SNe and AGB stars, 
stellar mass loss and SN feedback all take place on the 
timescale dictated by stellar evolution considerations. It is 
however important to note that published nucleosynthesis 
yields and stellar lifetimes vary significantly amongst au- 
thors. A detail e d discu ssion of this topic may be found in 
W iersma et al.1 d2009bl). The stellar population s ynthesis 
model PEGASE2 (|Fioc fc Rocca-Volmeran"gej|l997l ) is used 
to calculate the luminosity of the simulated galaxies. 

We now discuss the technical implementation of impor- 
tant "subgrid" processes included in our simulations. 



2.2 Supernova energy and mass loss 

In our simulations each star particle represents a simple stel- 
lar population (SSP), specified by its IMF, age, and metallic- 
ity. All stars more massive than 8Mq are assumed to explode 
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as SNe II. As the simulation proceeds, the supernova energy, 
mass loss, and associated metals produced by a star parti- 
cle are spread out over nei ghbouring parti c les. W e use the 
lifetime and yields given bv lPortinari et al l (|l998l ). While it 
would be most natural to distribute these quantities at each 
of the timesteps used to compute the dynamics, this proves 
to be prohibitively expensive computationally (particularly 
for SNe la). We adopt instead a more econom ical stochastic 
scheme, introduced bv lOkamoto et aD l|2008bl ). in which the 
distribution operation takes place within a smaller number 
of discrete time intervals during the lifetime of an individual 
star particle. 

From the time a star particle is born until it reaches 
the age, tsM Q , corresponding to the lifetime of an 8Mq star, 
(~ 40 Myr for solar metallicity) , the distribution of mass, 
metals and supernova energy to neighbouring particles oc- 
curs in a series of discrete events set by a timestep, Atu, of 
magnitude a fiftieth of t$M e - After an age tsM e there are no 
more SNe II, and a longer controlling timestep is introduced: 
Ati a = 100 Myrs. This longer step is chosen specifically to 
follow SNe la, but it is also sufficiently small to allow the 
distribution of elements produced by AGB stars, as well as 
stellar mass loss from the SSP as a whole. In general, the 
dynamical timestep, At, is much shorter than either Atu or 
Ati a . As a precaution, however, for star particles with an 
age less than tgM e , the dynamical time is explicitly limited 
to a maximum of half of Atu- 

For a star particle with an age, t < t$M Q , we define 
the conditional probability that the star particle undergoes 
a discrete event resulting in the distribution of energy, mass, 
and metals during a particular dynamical timestep, At, as 
follows: 



f t t+At m(t')dt' 



Pn = 



/; o 0+AtlI r„(t')dt'-// o ni(t')dt' 
j; +At ru(t')dt> 



t +At„ 



rn(t')dt' 



(1) 



where m(t) is the SN II rate for the SSP of age t and to 
is the age when the previous SN II timestep finished. The 
SN II rate, m(i), can be directly computed from the IMF, 
$(m) = d7V/dm, and the stellar lifetime, t(m) as 



. . , . dm(t) . . , 

ra[t) = -<P(m) for t(m u ) < t < t S M Q , 

where the IMF is normalised as 

m$(m)dm = 1, 



(2) 



(3) 



and m u and m; are the upper and lower limits of the IMF 
respectively. (We assume m u = 100 M Q and mj = 0.1 M§) 
The total number of SNe II produced by an SSP of 1 Mq is 
1.18x 10 -2 . We generate a uniform random number between 
zero and one at every At. When pn exceeds this number, the 
energy, mass, and metals expelled by the SSP during the 
interval to to to + Atu are distributed over the neighbouring 
gas particles and the probability is set to zero until the end 
of this SN II timestejB 

The energy, metals and mass produced by a particle are 



For example, the energy released by an SSP of IMq during 



distributed amongst its 40 nearest neighbouring gas particles 
at each of these discrete events. It is important that each 
quantity that is being redistributed be accurately conserved. 
To ensure that this happens, the fraction of each quantity 
received by the i-th gas particle neighbouring the star is 
given by rrn/^,. rrij, where the subscript j runs over all 40 
neighbours. We find that the results are relatively insensitive 
to the exact number of neighbours, with 32-128 particles 
producing similar behaviour. 

The same procedure applies after an age tsM e , except 
that the timestep is now Ati a . We calc u late t he SN la rate 
using the scheme of Grcggio & Renzini (1983) with the up- 
dated parameters of lPortinari et al.l (|1998l ). We take chem- 
ical yields and st e llar m ass loss from intermediate mass 
stars from Marigo (2001) . Further details may be found in 
lOkamoto et al.l (120051 ). iNagashima et al.l (120051 ) and refer- 
ences therein. 

In summary, our method for treating feedback and mass 
loss allows us to perform the expensive neighbour search 
only once per star particle per SN timestep and results in 
the correct distribution of SNe in time, even with large SN 
timesteps, when averaged over many particles. Note that 
the method guarantees a single distribution event per star 
particle per SN timestep. 



2.3 Galactic winds 

In our subgrid model of feedback, SN explosions give rise to a 
wind by imparting kinetic energy to nearby gas particles. We 
assume that all of the energy released from SNe is potentially 
available to power the kinetic energy of the wind. Hence in 
terms of feedback efficiency, there is no difference among our 
models. 

In order to prevent gas particles outside star-forming 
regions from being launched as windifl, we assume that only 
star- forming gas particles (tih > nn,th = 0.1 cm -3 , where 
nji,th is the threshold density above which star formation 
can occur) are eligible to become wind particles. If low den- 
sity gas particles (nu < fiH.th) receive feedback energy, we 
simply add it to their internal energy. 

During any given timestep, an eligible gas particle may 
receive supernova energy, AQ, from one or more neighbour- 
ing star particles. If this happens, then the particle is se- 
lected to become a wind particle during that timestep with 
a probability: 



Pw = 



AQ 



(4) 



where ttisph is the mass of the gas particle and n w is the 
initial wind speed. The value of t> w is discussed later in this 
subsection. Having decided whether a particle is to become 
a wind particle or not, we set its value, AQ = 0. 

This procedure ensures that the kinetic energy is in- 
jected locally and that the generation of winds is consistent 
with the timed release of SN energy and metals. This is 



Atu is Ssn Jt® +Atu m(t')dt', where Esn is the energy released 
by each supernova which, throughout this paper, we take to be 
E S n = 10 51 erg 

4 This can occur, for example, if a SN la is produced by a star 
in a low density gas environment. 
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an important difference f rom th e widely used wind model 
of ISpringel fc Hernquistl (|2003al ) in which the wind parti- 
cles are selected stochastically from all the star-forming gas 
(i.e. with nu > nn,th) in the simulation and are there- 
fore not local to young star particle^]. A different im- 
plementation of locally gener ated winds can be found in 
iDalla Vecchia fc Schavel feoosl ). 

It is possible for p w to exceed unity. We deal with this 
situation using an iterative procedure. Firstly, the particle is 
added to the wind and we record both its label and the value 
of its 'excess' energy, AQ — |mspHVw- Having considered all 
eligible gas particles to decide which become wind particles, 
we then revisit those with excess energy and redistribute 
this energy to neighbouring non-wind gas particles. We then 
consider those particles whose density makes them eligible 
to be selected as a wind particle, recompute p w including 
the reassigned excess energy from AQ in Eqn. [4] and decide 
whether they are converted into wind particles or not in the 
usual way. If there are still values of p w which exceed unity, a 
further iteration is begun, and these iterations continue until 
there is no more excess energy to redistribute. In practise no 
more than two iterations are ever required. 

When a particle is selected as a wind particle, its veloc- 
ity, v, is incremented as: 



v + v w n. 



(•») 



The unit vector n is chosen at random to be either par- 
allel or anti-parallel to the vector (v — v) x a grav , where 
v ideally would be chosen to be the velocity of the cen- 
tre of mass of the galaxy. Since it is difficult to determine 
this velocity as the simulation proceeds, we instead select 
the average velocity of the neighbouring 40 dark matter 
particles as a proxy. The gravitational acceleration, a grav , 
on the other hand, includes contributions from matter ev- 
erywhere in the simulation. In practise, for the central re- 
gions of haloes and subhaloes, both the magnitude and di- 
rection of a grav are largely determined by the local mass 
distribution. Adding a velocity according to Eqn. [5] leads 
to wind particles being ejected preferentially along the rota- 
tion axis of a spinning object, thus generating an 'axial wind' 
|Springel fc Hernqui st 2003a). Note that while the momen- 
tum is not strictly balanced when a single particle is placed 
in the wind, the momentum and angular momentum are sta- 
tistically conserved because of the random orientation of the 
momentum kick along direction n. 

Since the physical mechanisms that drive the outflows 
in galaxies are not well understood and the mass in winds 
is hardly constrained by observations, we adopt two phe- 
nomenological models for the properties of the winds (both 
of which assume that all of the supernova energy is available 
to power the outflows). The first model (which we label W 
for 'variable wind') assumes that the initial wind speed, v w , 
is proportional to the one-dimensional velocity dispersion, 
a, determined from the 40 neighbouring dark matter par- 
ticles of the gas particle in question. That is, in the 'vw ' 
model, v w tx a, as suggested by recent data l|Martinll2005h . 



5 Note however that in ISpringel fc Hernquistl 1 1200331 1 the gas par- 
ticle metallicity is updated according to its star formation rate; 
hence the non-local winds are consistent with the chemical en- 
richment scheme. 



Here, we are taking the local velocity dispersion as a proxy 
for the circular velocity. We have checked that the local ve- 
locity dispersion measured at the position of a star-forming 
gas particle is strongly correlated with the maximum of the 
circular velocity profile of its host (sub)halo, Vmax, and that 
the relation between them (Vinax — 1.45er) does not evolve 
with redshift. 

In order to characterise the 'vw' models, we discuss 
wind properties by assuming the instantaneous recycling 
approximation and considering only SNe II as the energy 
source for the winds. While we do not make these assump- 
tions in the simulations, considering this simplest of cases 
is instructive. With these assumptions, the energy ejection 
rate, Esn, is proportional to the star formation rate, M,. 
Therefore, the mass flux in the winds in the 'vw' models, 
M w , can be written as: 



M„ = 2—- tx — — 

V <T 



because the wind speed, u w , is proportional to the local dark 
matter velocity dispersion, a. From these considerations, the 
wind mass- loading, 7? w , can be expressed as: 



(7) 



where ao is a parameter defined by the IMF and the feedback 
efficiency. We give the values of ao in our models in Table [2] 
As mentioned earlier, we adopt the Chabrier IMF and make 
use of all the SN energy to power the kinetic energy of the 
winds. 

For the second wind model (which we label 'cw' for 
'constant wind'), we assume a constant initial wind speed, 
which implies a constant mass-loading. This type of wind 
model has been widely used in hydrodynamic simulations 
of galaxy formation (e.g. Springel & Hernquist 2003a; 



i Nagamine. Springel. fc Hernquistl |2004| ; lOkamoto et al.l 
l2008bl ). 

In both the 'vw' and 'cw' models, we assume that newly 
launched wind particles are decoupled from hydrodynamical 
interactions for a brief period of time in order to enable 
the winds to emerge from a location close to the surface of 
their star-forming region with the specified wind velocity and 
mass-loading. This wind decoupling h as been widely used in 
simulat ions since it was introduced bv lSpringel fc Hernquistl 
(l2003ah 



Dalla Vecchia fc Schavel (|2008l ) compared coupled (i.e. 



not decoupled) and decoupled winds in detail and found that 
coupled winds blow bubbles and drive turbulence or create 
channels in the gas disc, while decoupled winds remove fuel 
for star formation without disturbing the gas disc. More 
importantly, for a given set of (input) wind speed and mass- 
loading, coupled winds are more efficient at suppressing star 
formation in low mass galaxies because they drag gas along; 
consequently the mass-loading when they leave the star- 
forming region (i.e. the effective mass-loading) is much larger 
than the input one. On the other hand, coupled winds are 
less efficient in high mass galaxies because they suffer large 
energy losses due to drag in the high-pressure ISM and the 
decelerated winds cannot escape fr om the galaxy (or often 
even from the star-forming region) (|Schave et alj||201Ch . 

Since we wish to investigate the effects of two differ- 
ent feedback models on the population of satellite galaxies, 
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in which wind speed and mass-loading depend differently on 
the galaxy properties, we employ decoupled winds in this pa- 
per in order to specify the effective wind velocity and mass- 
loading as inputs. Full hydrodynamical interactions are en- 
abled once a wind particle leaves the star-forming region 
(Vih < O.lriH.th) or after the time, 10 kpc/w w , has elapsed, 
whichever occurs earliei0. 



2.4 Models 

For the variable wind models, we show results for u w = 
5cr and v„ — 4a; we refer to these as 'vw5<r' and 'vw4cr', 
respectively. Although we also calculated a model with v v — 
3cr, we do not show the results here because this wind speed 
is insufficient to allow gas particles to escape from their host 
halo and thus star forma tion is hardly suppress ed compared 
to the no feedback case (|Okamoto et al.ll2008bl ). 

For the constant wind models, we show results for 
u w = 700 and 600 km s _1 ; we label them 'cw700' and 
'cw600', respecti vely. We employ a slight ly faster wind speed 
than was used in lOkamoto et al.l l|2008bl ) (v w = 500 km s _1 ) 
because we find that this wind speed is too low to suppress 
the early starbursts that would turn our g alaxies into bulge- 
dominated objects (|Okamoto et alj|2005h . We believe that 
this is due to the fact that the baryon fraction in our cur- 
rent model s is 35% higher than in the cosmological model 
adopted bv lOkamoto et all |2008b). The values of the wind 
mass-loading factor, r/ w , for these models are given in Ta- 
bled 

For the models introduced so far, we employ a mul- 
tiphase model for the star-forming g described in 

lOkamoto et~ai1 (|2008bh . in which a minimum pressure for 
star-forming gas, P m in oc p 1A , is imposed in order to stabilise 
gas discs against gravitational instability. Such a stiff equa- 
tio n of state for the star-formin g gas was first introduced 
bv lSpringel fc Hernquistl l|2003af ) in an attempt to allow for 
the fact that the interstellar medium (ISM) is supported 
by a hot gas phas e which is unresolved in the s imulations. 
On the other hand. lWada fe Norman] i|200ll . l2007l ) suggested 
that the multi-phase ISM is supported by turbulent motions 
induced by the self-gravity o f the g as and galactic rotation. 
In any case. lRobertson et al.l ((2004) argued that a stiff equa- 
tion of state for the star-forming gas is required to form a 
disc. Our adoption of such an equation of state is intended 
to allow for these unresolved and ill-understood processes. 

Artificial fragmentation of a disk may occur when 
simulat ions do not properly resolve the loca l Jeans 
length l|Truelove et all 1 19971 ; iBate fc Burkertlll997f ). In or- 
der to prevent this from happening, several authors have 
adopted a similar strategy to us by imposing a minimum 



pressure for high- 



Saitoh et al.l 



Schave et al 



d ensity gas (jRobertson fc Kravtsovl 



2008; Co vering et al.l 120091 ; ICrain et all 



2010). 



2008 



2009; 



Schave et alJ(|201Ch found that the global star forma- 



tion rate is almost completely independent of both the equa- 
tion of state of the ISM and the star formation law because 
the star formation is highly self-regulated by the winds. 
We therefore perform simulations without this multiphase 



B In our simulations, the first condition (njj < 0.1riH,th)i is usu- 
ally satisfied first. 



Table 2. Wind models and their properties. ti w and r; w denote 
the initial wind speed and the wind mass-loading (see text). 



Model 



v w (km s 1 ) 



'/w 



Multiphase 



vw5(T aingle 
vw5cj 

vw4cr 

cw700 
cw600 



5(7 
5(7 
4(7 

700 
600 



217km 



v 271km s- 1 , 

2.4 
3.3 



No 

Yes 

Yes 
Yes 
Yes 
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Figure 1. Distribution of gas particles in the density- 
temperature phase diagram for the Aq-A halo at z = 0. The 
left and right panels show the 'vw5(7 s ; n gi e ' and 'vw5(7' models, 
respectively. The density above which star formation is enabled 
is indicated by the vertical dotted line. The ISM is well enriched 
by z = 0. As a result there is low-temperature (< 10 4 K) star- 
forming gas in the 'vw5(7 s ; ng i e ' model which has been cooled by 
metal-line cooling. We show the effective temperature for the mul- 
tiphase star-forming gas in the l vw5(7' model. Most of the star- 



forming gas lies on the line specified by P n 



oc p 



model for comparison. In this case, the star-forming gas can 
cool down to the equilibrium temperature set by the bal- 
ance between photoheating and radiative cooling. Note that 
metal enriched gas can cool far below 10 4 K b ecause we in- 
clude metal- line cooling ijWiersma et al.ll2009ah . We use this 
single-phase gas in a 'vw5cr' wind simulation and label it 
'vw5(7 sing i c '. 

In Fig. Q]we compare the gas distribution in the density- 
temperature phase diagrams for the 'vw5(7 s i ng i e ' and 'vw5cr' 
simulations in the Aq-A halo at z = 0. Whereas the star- 
forming gas in the 'vw5cr s i ng ie' model cools down to ~ 100K 
because of metal-line cooling and an increasing cooling rate 
with density, the gas in the 'vw5a' model lies on the imposed 
minimum temperature. The equation of state of the star- 
forming gas in the 'vw5a' model is hence much stiffer than 
that in the 'vw5cr a i ng i e ' model. For both ISM models, the star 
formation efficiency is normalised to reproduce the observed 
relation between the sur face star formati on rate density and 
the surface gas density l|Kennicuttlll99&f ). 



3 RESULTS 

We begin this section by presenting an overview of the 
properties of the central galaxies. While our main focus 
is on the satellite galaxy population, the bulk properties 
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Table 1. Numerical parameters relating to our four sets of initial conditions. The virial mass, M v ; r , is denned as t he mass within a 
sphere whose mean density is equal to the virial density at redshift zero computed from the spherical collapse model l|Eke et al.| [l996). 
The gravitational softening lengths, e, are kept fixed in comoving coordinates for z > 3; thereafter they are frozen in physical units and 
to values presented in the table. 



Halo 


M vir (fc- 1 ^) 


m D M (h 1 M Q ) 


msPH (k -1 M Q ) 


e (h- 1 kpc) 


Aq-A 


1.4 x 10 12 


1.9 x 10 6 


4.1 x 10 5 


0.425 


Aq-C 


1.2 x 10 12 


1.5 x 10 6 


3.4 x 10 5 


0.425 


Aq-D 


1.3 x 10 12 


1.9 x 10 6 


4.1 x 10 5 


0.425 


Aq-D-HR 


1.3 x 10 12 


1.6 x 10 5 


3.5 x 10 4 


0.175 



of the central galaxies, including their morphology, are a 
useful aid in understanding the nature of our wind models 
which, in turn, are the key to understanding the properties 
of the satellites dMaller fc Dekelll20o3 ; lOkamoto etai1l2005l ; 
IZavala et al.ll2008l ). A more detailed analysis of the central 
galaxies will be presented in a forthcoming paper. 



3.1 The central galaxies 

A visual impression of the structure of our simulated galax- 
ies may be gained from Fig. [2] which shows edge-on views 
of the galaxies at z — 0. To make these images, we define 
a l z' axis for each central galaxy as the direction of the 
angular momentum of all stars within 0.05J? v ir of the cen- 
tre at z = 0. We take 7? v ir to be the radius of a sphere 
that has the virial overdens ity given by the spherical col- 
lapse model (|Eke et al.lll996h . We define two mutually per- 
pendicular ' and l y\ at random for each simulation, 
perpendicular to the z-axis. Fig.[2]shows that striking differ- 
ences in the morphologies of the galaxies arise in each of the 
haloes Aq-A, Aq-C and Aq-D, when the wind presc r iption 
is varied. This supports the claim of iOkamoto et all (|2005l ) 
that the morphology of simulated galaxies is very sensitive 
to the details of the physics implemented in the simulation. 
For example, for halo Aq-C a mere change of ~ 15% in 
the wind speed from 700 km s _1 to 600 km s _1 results in a 
change from a disc-dominated galaxy to a spheroidal galaxy. 
Nonetheless, it is encouraging that some models do success- 
fully form disc-dominated galaxies. 

The morphology of a galaxy is largely dictated by the 
angular momentum distribution of the stars. To quantify 
this, we analyse the distribution of the orbital circularity, 
denned below, of stars in the central galaxies. We exclude 
stars belonging to the satellites (which we define precisely 
in the next subsection). For each star belonging to a central 
galaxy we compute, J z , the component of specific angular 
momentum parallel to the z-axis defined earlier. We then 
compute the specific angular momentum, J C (E), of a pro- 
grade circular orbit with the same binding energy as the 
particle. To evaluate J C (E), we first compute the circular 
velocity profile of the host halo, v c (r). The specific binding 
energy of a particle on a circular orbit at radius r is then 
given by 



GM(< r) 1 . . 2 
E(r) = i L + - Vc ( r f. 





OK 



cw600 



Figure 2. Edge-on views of the central galaxies at 2 = 0. The 
brightness is coded by the B-band surface brightness within a 
50/i kpc box centred on each galaxy. The viewing angles are 
determined using the angular momentum of stars within 5% of the 
virial radii of the host haloes. The Aq-A, Aq-C, and Aq-D haloes 
are shown from left to right, and models, 'vw5<7 s i ng i c ', 'vw5ct', 
l vw4o-', l cw700', and l cw600' are shown from top to bottom. 



(8) 



Similarly, the specific angular momentum of a particle at r 



J c (r) = rv c (r). 



(9) 
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Figure 3. The mass-weighted normalised distribution of orbital 
circularity, J Z /J C (E), of all stars at z = lying within i? v ir of the 
centre, excluding stars identified as members of satellite galaxies. 
The Aq-A, Aq-C, and Aq-D haloes are shown from left to right 
and the 'vw5<r', 'vw4cr', and 'cw' models are shown from top to 
bottom. A peak at J z /jc(E) ~ 1 indicates the existence of a disc 
component. The '™5o- s i ng i e ' and 'cw700' models of the Aq-A halo 
and the 'vw5cr' model of the Aq-C halo exhibit a retrograde disc 
component as well as a prograde one. 



We then calculate the specific energy of the i-th star particle 
assuming a spherically symmetric halo: 



Ei 



GM(< n) , 1 2 
n 2 



(10) 



We identify the radius, r, of the circular orbit of spe- 
cific binding energy, E(r) = Ei using a Newton-Raphson 
method. This then gives the value of J c (Ej). The ratio 
Jz/Jc(E) defines the orbital circularity (|Abadi et al.l l2003b; 
lOkamoto et al.ll2005T ). 

Fig. [3] shows the normalised distribution of the orbital 
circularity of stars in the central galaxy within R V1I . A cold 
disc component has J Z /J C (E) ~ 1 and such a component is 
evident in most of the galaxies except in the 'cw600' model 
of the Aq-C halo. The 'cw600' model produces the most 
disc-dominated galaxy in the Aq-A halo, but the least disc- 
dominated galaxy in the Aq-C halo. For a fixed wind model, 
halo-to-halo differences in the final galaxy morphology re- 
flect differences in the formation and merging histories of 
the haloes (even though most of them have similar quiet 
histories at early times). Galaxy morphology is particularly 



sensitive to formation history in the 'cw' models because of 
the introduction of a critical velocity scale for the windfl 
In the W case, the dependency on halo formation history 
is much weaker. 

In Aq-A and Aq-D, the most pronounced discs among 
the 'vw' models form in the 'vw5cr' model, but in Aq-C, it is 
the 'vw4a' which produces the largest disc. The 'vw5a- s i n gie' 
galaxies are less disc-dominated than their 'vw5<r' coun- 
terparts. That confirms an earlier claim that the stiffer 
equat ion of state of the sta r-forming gas helps disc forma- 
tion (|Robertson et all 120041 ). The effect is, however, not al- 
ways evident (see Aq-C). Interestingly, the 'vw5o- s ingic' and 
'cw700' models of the Aq-A halo and the 'vw5er' model of 
the Aq-C halo exhibit counter-rotating or retrograde discs in 
addition to the prograde discs. The origin of these counter- 
rotating discs will be discussed in a forthcoming paper. 

Our 'vw' galaxies are more disc-d ominated than those 
simulated by IScannapieco et al.l |2009) in the same Aquar- 
ius haloes, using very similar initial conditions but a differ- 
ent galaxy formation code. Since galaxy morphology is very 
sensitive to the details of feedback, this discrepancy is not 
surprising and probably reflects the fact that the strong feed- 
back in the 'vw' models is able to suppress star formation in 
small haloes more effic iently than that in the ir simulations. 

As emphasised bv lOkamoto et all (|2005l ). the star for- 
mation history of a galaxy is correlated with its morphology. 
Star formation histories for the central galaxies in our sim- 
ulations are shown in Fig. [4] The 'vw5<r' and 'vw5o" s ingie' 
models are very similar, indicating that the star formation 
history is even less sensitive to the adopted ISM model than 
the morpholog y of the galaxy. T his is consistent with the 
conclusions of lSchave et al.l (|2010T ) who find that the global 
star formation is highly self-regulated by feedback and is al- 
most independent of the choices of equation of state for the 
ISM and star formation law. 

The 'vw5cr' and 'vw4a' models are also similar except 
that the 'vw4<r' model for the Aq-D halo has significantly 
higher star formation rate at low redshift. It is interesting 
that the 'cw600' model of the Aq-C halo has a relatively 
high star formation rate at low redshift, a feature which 
was c laimed to be a sign of disc formation ( Ok amoto et al.l 
2005), although in this case this galaxy is dominated by 
a spheroidal component (see Fig. and [3]). Note that our 
sample of Aquarius haloes all have qui et merger histories 
at low redshift. As recently stressed by IScannapieco et al.l 
(2009), the origin of galaxy morphology is more complicated 
than we generally thought and the star formation history is 
only a very crude indicator of the final morphology. 



3.2 The satellite galaxies 

We now turn to an investigation of how different feedback 
models affect the properties of satelli tes. We identify sub - 
haloes using the SUBFIND algorithm |Springel et alj|200ll ) 
and restrict attention to subhaloes of at least 32 particles. In 
counting the number of particles in a subhalo, we include not 
only dark matter but also baryonic particles. Therefore it is, 



7 This is not true f or coupled winds (see 
Dalla Vecchia & Schaveiliooih . 
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Figure 4. The star formation histories of the central galaxies. 
Only stars within 0.1i? v i r of the central galaxy which are not 
members of a satellite galaxies are shown. The panels and lines 
are the same as in Fig. [3] The 'vw5a- s inglc' and 'vw5cr' models 
show almost identical star formation histories, implying that the 
modelling of the ISM has little effect on the star formation his- 
tories and that it is the wind model which is the most important 
ingredient in determining the morphology of the galaxy. 



in principle, possible to identify purely baryonic gravitation- 
ally bound objects. We define satellite galaxies as subhaloes 
that contain at least 10 star particles. In Sec. 3.3 we check 
the robustness of this definition and carry out a convergence 
test. 



3.2.1 The satellite luminosity junction 

In Fig. [5] we plot the satellite luminosity functions in our 
simulations and compare them with observations of the 
Local Group. The observed luminosity function of bright 
[My < —11) satellites lying within 280 kpc of either the 
Milky Way or M31 is construct ed from the d ata compiled by 
Metz. Kroupa. fc Jerienl (120071 ). The data of lKoposov et all 
20081 . which includes the incompleteness-corrected SDSS ul- 
trafaint satellites within 280 kpc from the centre of the 
Milky Way, are represented by a power-law. The host 
;alaxy haloes, de fined by the friends-of-friends algorithm 
Davis et alii 19851 ). extend beyond 280 kpc, but, for consis- 
tency with the observations, in constructing the luminosity 
functions we consider only the satellites within 280 kpc of 
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Figure 5. Satellite luminosity functions. Only satellites within 
280 kpc from the central galaxies are counted. The panels are the 
same as in Fig. \3\ The solid lines show the results from the sim- 
ulations. The dashed line shows the luminosity function for the 
bright (My < -11) Milky Way and M31 satellites, similarly tak- 
ing only those within 280 k pc from eithe r galaxy, constructed from 
the data compiled bv lMetz et al.l ||2007| ). T he dot-dashed l i ne is a 
power-law fit to the luminosity function of Koposov et al. | i2008h 
which includes newly discovered SDSS ultra-faint satellites, again 
restricted to include only satellites within 280 kpc of the Milky 
Way. 



the centre. For other analyses below we will make use of 
larger samples of satellites. 

For a given wind model, the satellite luminosity func- 
tion varies significantly from halo to halo, refl ecting the vari- 
ation in the number of m assive subhaloes (|Springel et al.l 
2008; Ishivama et alJ l2009h . This alerts us to the danger 
of assuming that the Milky Way is a typical galaxy. The 
'vw5<T s ingie' and 'vw5cr' models again yield almost identical 
results. Interestingly, the 'cw700' and 'cw600' models also 
produce results almost indistinguishable from each other. 
The 'vw5cr and 'vw4<t' models are similar except that the 
brightest satellite in a 'vw4rr' model is slightly brighter than 
its counterpart in the corresponding 'vw5a' model reflecting 
the lower wind speed in the former model. 

Our results are co nsistent with the co nclusion of the 
semi-analytic study by iBenson et al. (2002) that it is diffi- 
cult to form satellite galaxies as bright as the LMC in Milky 
Way-sized haloes. Although in the 'vw4<r' model for the 'Aq- 
D' halo there is a nearby galaxy as bright as M33, it does 
not qualify as a satellite galaxy according to our strict deli- 
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nition which requires satellites to be within a radius of 280 
kpc from the central galaxy. 

Away from the bright-end, the 'vw' models (in which 
the wind mass-loading increases towards smaller subhaloes) 
do reproduce the observed satellite luminosity function rea- 
sonably well. By contrast, both of the 'cw' models have 
overly steep luminosity functions with too many faint satel- 
lites. The constant mass-loading wind model is not able to 
suppress star formation in these faint satellites sufficiently. 
We expect that if we had adopted the mom e ntum-drive n 
wind model favoured bv lOppenheimer fc Dave] |2006, 2008), 
in which the wind mass-loading, n v , is proportional to cr _1 , 
we would have obtained satellite luminosity functions whose 
faint-end slopes are steeper than those in our 'vw' models 
(rjy, oc &~ 2 ). On the other hand, it seems likely that if we 
had not assumed that the winds are decoupled for a short 
period time, the faint end of the resulting luminosity func- 
tions in the 'cw' cases woul d have been shallower than in 
our decoupled 'cw ' models (|Dalla Vecchia fc Schavd [20081 : 
ISchave et afll201Ch . Within any given halo, the 'vw' models 
always match the satellite luminosity function better than 
the 'cw' models. From this point of view there is little to 
choose between the 'vw5<r' and 'vw4cr' versions. 

The robustness of our results depends on how well the 
low-mass end of the subhalo mass function is resolved in 
the simulations. In order to identify the smallest resolved 
subhalo, we first compute the cumulative subhalo abundance 
as a function of maximum subhalo circular velocity, 7V(> 
Vmax). We then fit this cumulative velocity function with a 
power- law, N(> Vmax) oc Vm ax ; subhalo velocity functions 
in N-body simula tions are well represented by a power-law 
with a - -3 (e.g. lDiemand et al.ll2007l ; ISpringel et alj|2008t ) 

We define the smallest resolved circular velocity as the 
Vmax value at which the velocity function first deviates by 
10% from the power-law fit. We have computed this for the 
'vw5<r' and 'cw700' models of the Aq-D halo, which is the 
halo we use for convergence tests later. For the 'vw5er' model 
we find a value of 12.9 km s _1 , which is consistent with 
the V m ax at which the subhalo velocity function begins to 
deviate from that of the high-resolution counterpart (Aq- 
D-HR). The smallest V max of the subhaloes which contain 
at least 10 star particles (i.e. the satellites) is 22.8 km s _1 . 
We conclude that the subhaloes that host satellites in the 
'vw' models are well resolved. For the 'cw700' model, the 
minimum resolved circular velocity is 11.2 km s _1 , which 
is consistent with the value for the vw5<r model. However, 
now the smallest V m ax of the satellites is 13.2 km s" 1 , which 
is close to the resolution limit. Thus, in the 'cw' models, 
the flattening seen at the very faint end of the luminosity 
functions could be affected by numerical limitations. 



3.2.2 Metallicity and abundance patterns 

We now investigate the metallicity of satellites which, 
in principle, is also a sensitive diagno s tic of the feed- 
back physics (e.g. lArimoto fc Yoshn]|l987l ; iFinlator fc Dave] 
2008). In Fig. IS] we plot the iron abundance as a function of 
satellite luminosity. (To compare our simulations with the 
data, we adopt the values for the solar abundances given 
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Figure 6. The iron abundance relative to solar for satellite galax- 
ies, as a function of the absolute V-band magnitude. The panels 
are the same as in Fig. [3] The abund ances for the Local Group 
satellites are taken from lMateol l|l998t ). wi th the exception of the 
Magellanic Clouds which was taken from lffill et~aH (Il995ll . Dif- 
ferent symbols correspond to different wind models as indicated 
in the legend. The Local Group data are shown by plus signs. 

by iGrevesse fc Sauvall (|l998l f § l') All of the 'vw' models re- 
produce the observed luminosity-metallicity relation for the 
Local Group satellites quite well. On the other hand, the iron 
abundance of the 'cw' satellites is almost a constant func- 
tion of the satellite luminosity, reflecting the constant wind 
mass-loading factor. These results strongly support mod- 
els in which the wind mass-loading scales as ry w oc a~ 2 . A 
momentum-driven wind would most likely produce a rela- 
tion intermediate between those of the 'vw' and 'cw' mod- 
els, thus too shallow to match the data, while a coupled 
wind model would like produce a steeper relation than the 
corresponding 'cw' and 'vw' models. 

The alpha element-to-iron abundance ratio is often 
used as a measure of the s tar formation timescale (e.g. 
iNaeashima fc Okamotd l2006i. Since SNe la hardly produce 
alpha elements, [O/Fe] decreases with time if a star for- 
mation episode lasts longer than ~ lGyr, which is a typ- 
ical SNe la timescale. For an SSP of solar metallicity, it 
takes ~ 0.9 Gyr for SNe la to produce as much iron as 
that produced by SNe ifl In Fig. [7] we show the oxygen- 

8 Assuming the more re cent values for the solar abundances given 
bv lAsplund et alj [|2005h increases [Fe/H] in the simulations only 
by ~ 0.16. 

9 Assuming solar metallicity, the first SN la goes off at t ~ 75 
Myr in our model. 
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Figure 7. The oxygen-to-iron abundance ratio of satellite galax- 
ies as a function of the iron abundance. The panels correspond 
with those in Fig. \E\ In addition to data from the sources given 
in Fig. we plot the median values for Local Group dwarf 
spheroidals taken fr om IShetrone, Cote, fc Sargentl J200lh and 
IShetrone et al.l i2003l) . For Mateo's data, the error bars represent 
1(7 error, while for the Magellanic Clouds and dwarf spheroidals 
they show the larger of the 25-75 percentile interval or the obser- 
vational 1<t error. 



to-iron abundance ratio of satellite galaxie^j as a function 
of their iron abundance. Most of the satellites in the W 
models have [O/Fe] ~ 0, which is slightly lower than the 
observed values. This might imply that the star formation 
timescales in the simulated dwarf satellites are too long, or 
that the model emplo yed for the SN la is incorrect (e.g. 
Kobayashi ct al. I ll998l ; o ur SN la rate was calcula ted using 
the scheme described bv lGreggio fc Renzinilll983l. with pa- 
rameters updated according to Portinari et aljl998h. or tha t 
the assumed IMF is inaccurate fe.g. lNagashima et alJ r2005). 
Both the SN la rat e and the nucleosynthes is yields are highly 
uncertain (see e.g. IWiersma et al.ll2009bl ). Exploring these 
issues is, however, beyond the scope of this paper and we 
defer this to future work. 

The 'cw' models reproduce the observed abundance 
patterns in the Milky Way satellites reasonably well, even 



10 Again, using more recent solar abundance determinations de- 
creases [O^Fe] in the simulations by 0.17. However, the data in 



iMateol dl99Sl) assume 12 + log(Q/H) w = 8.93, whic h is consistent 
with the values given bv lGrevesse &: Sauvall lll998Tl 



though they fail to match their iron abundance (see Fig. [6j) . 
The reasons behind the different patterns in the W and 
'cw' are instructive. In the W models, the wind mass- 
loading is larger for smaller galaxies resulting in more metal 
poor dwarfs (see Fig. [6j ; this, in turn, suppresses star forma- 
tion. However, the wind speed, t> w ~ 4-5cr, allows ex pelled 
gas to fall back later fsee lOppenheimer fc D ave 2008), thus 
lengthening the effective star formation timescale. In the 
'cw' models, there is a critical velocity at which the wind 
speed is equal to the escape velocity from the halo, v csc . For 
l?w Vcsc, the gas blown out of a halo never comes back 
and therefore the winds shorten the effective star forma- 
tion timescale. As a result, galaxies with small F max become 
alpha-enhanced. When the wind speed is of the order of the 
escape velocity, the expelled gas is recycled as a fountain. 
Hence, when « w < v aS c, the effective star formation timescale 
is longer than when v w 3> v CS c and the oxygen-to-iron ratio 
begins to decrease relative to the value set purely by SNe II. 



3.2.3 Star formation histories 

The star formation histories of randomly selected satellites 
in two models, 'vw5<j' and 'cw700' in the Aq-D halo, are 
shown in panels (a) and (b) of Fig. [8] respectively. (We have 
checked that the other models of each type are similar to 
these.) These plots allow us to investigate how reionization 
and infall into a larger halo affect star formation. The ar- 
rows indicate the 'accretion' redshift, i.e. the moment when 
the main progenitor of a satellite actually becomes a satel- 
lite galaxy by falling into a larger halo. The labels give the 
maximum of the circular velocity of the selected satellites, 
Vmax, in km s -1 , at z = 0. 

We find that reionization (which happens at z IC = 9) 
has little impact on the star formation histories of today's 
satellites. Even the smallest of t hem continue to ma k e star s 
past this epoch. As shown by Oka moto fc Frenkl (2009), 
however, this is not the case for haloes that are too small 
at the epoch of reionization to host satellites. Haloes which 
have Knax < 12km s _1 at the epoch of reionization are un- 
able to hold on to their gas which is evaporated by the UV 
photons and thus end up as dark subhaloes today. 

The figure shows that infall into a larger halo is im- 
portant although it do es not always immediately truncate 
star fo rmation (see also lFont et al . 20081. lOkamoto fc Frenkl 
(2009) demonstrate that infall tends to strip gas off the in- 
falling small object, eventually choking its star formation. 
A decrease in star formation soon after infall is apparent in 
some of the examples in the figure, such as the haloes with 
s c m " = 41.6, 39.2, 28.5, 23.8 km s" 1 shown in panel (a) for 
'vw5ct'. In a few cases, e.g. w™ ax = 40.7, 31.4, 25.3 km s _1 , 
star formation does cease completely upon infall. 

Similar remarks can be made regarding the satellites in 
the 'cw700' model. In this case, because the wind mass- 
loading is smaller than in the 'vw' models, subhaloes of 
smaller circular velocity can host satellites. The truncation 
of star formation at the time of accretion is more evident in 
satellites with such small circular velocities. As may be seen 
by comparing panel (a) and (b) of Fig. [8] at a given Knax, 
the star formation timescale is shorter in the 'cw700' model 
than in the 'vw5er' model. In addition, the star formation 
in the l vw5o"' model is more variable in time than in the 
'cw700' model. These differences arise from the fact that a 




Figure 8. (a): The star formation rate divided by M®, the stellar mass at z = 0, for a selection of satellite galaxies in the vw5<r' model 
of the Aq-D halo. The labels give the value of Vmax of each of the individual satellites, in km s _1 . The arrows indicate the redshift when 
the progenitors are accreted into a larger halo (i.e. when they become satellites). Note that we assume reionization occurs at l + z re = 10. 
(b): As (a) but for the l cw700' model of the Aq-D halo. 



wind speed of v w = 4 — 5a in the W models can cause 
wind particles to escape from a halo which, however, can 
be re-accreted later on. By contrast, in the constant wind 
speed model, re-accretion is inhibited whenever v csc <C u w 
while a wind is suppressed whenever v csc > u w . As a re- 
sult, the star formation timescale is shorter than in the 'vw' 
models whether v CS c <C v w or v csc > v„. These differences in 
star formation histories explain the differences in abundance 
patterns seen between the two types of models (Fig. 0). 



3.3 A high resolution simulation 

In this subsection, we select the most promising wind model 
and investigate its sensitivity to numerical resolution by res- 
imulating it with ten times higher mass resolution. We se- 
lected the 'vw5cr' model on the grounds that it provides 
the best overall match to the observed luminosity function 
and iron abundances of Local Group satellites (although 
the oxygen-to-iron ratios, which are sensitive to the adopted 
SN la rate, are systematically low). The run, Aq-D-HR, is 
detailed in Table [T] (thanks to its low concentration, the Aq- 
D halo is the least computationally expensive of the three). 
Since we only consider the 'vw5<t' model in this subsection, 
for simplicity, we refer to the high and low resolution sim- 
ulations as Aq-D-HR and Aq-D (or simply as the high and 
low resolution simulations), respectively. 

We first examine the correspondence between the posi- 
tions of individual satellites in the two simulations. To do 
this, we trace the dark matter particles of a particular sub- 
halo in one of the simulations back to the initial conditions 
and search for a matche d subhalo in the other simulation 
fsee lSpringel et ai1r20 08j). We include all the satellites lying 
within 280 kpc in the low resolution simulation and the satel- 
lites in the high resolution simulation which also lie within 
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400 
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Figure 9. The projected positions of the matched satellites in 
the two simulations with different resolution, Aq-D-HR and Aq- 
D. The matching was carried out for satellites lying within 280 
kpc from the centre that are brighter than My = —11.2, cor- 
responding to the luminosity of the faintest satellite within 280 
kpc in the lower resolution simulation. Matched pairs are joined 
by thin lines. Circles denote satellites in Aq-D-HR and stars in 
Aq-D; open symbols indicate satellites found outside the 280 kpc 
sphere which is marked by the dot-dashed line. 
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280 kpc and are brighter than Mv = —11.2, corresponding 
to the faintest luminosity of satellites within 280 kpc in the 
low resolution simulation. 

The result is shown in Fig. [5] where the matched 
pairs are linked by a line. The mean spatial offset between 
matched pairs is ~ 100 kpc, somewhat larger than the mean 
offset of ~ 54 kpc found between m atched pairs in Aq-A-4 
and Aq-A-lE3 (ISpringel et alj 120081 '). This difference is not 
due to small number statistics - we have only 12 pairs of 
satellites - since a similar offset is found for matched pairs 
of all subhaloes within 280 kpc, whether or not they host 
satellites. The mismatch is more likely due to slight orbital 
phase deviations - which translate into large spatial separa- 
tions - produced by differences in the potential of the two 
galaxies whose evolution is not identical, as we will show in 
Fig. 1101 Slight difference in the mass of the subhaloes due 
to baryonic processes may also contribute to this offset. In 
spite of the relatively large spatial offsets, the photometric 
and chemical properties of the satellite populations in the 
two simulations are remarkably similar, as we will show later 
in this subsection. 

In Fig. 1101 we present the results of the convergence 
tests. We first compare the star formation histories in the 
central galaxies in two simulations in the upper-left panel. 
Although the shape of the two curves is very similar, 
the higher resolution simulation has a slightly lower mean 
star formation rate. This likely reflects enhanced mass loss 
through winds from small haloes not resolved in the lower 
resolution simulation. At z = 0, the stellar mass in the cen- 
tral galaxy is 18% lower in the Aq-D-HR simulation than in 
the Aq-D simulation. 

The upper-right panel shows the satellite luminosity 
functions obtained in the two simulations. They agree very 
well. The slight differences between them arise mainly from 
the requirement that the samples should contain only satel- 
lites found within 280 kpc of the central galaxy. Two of 
the satellites in the high resolution simulation have counter- 
parts in the low resolution simulation which lie outside the 
280 kpc radius. One of these can be clearly seen in Fig. [9] 
the other one is projected inside the circle. When we include 
these two satellites, the luminosity functions in the low and 
high resolution simulations agree even better, down to the 
faintest satellite luminosity in the Aq-D (Mv — —11). This 
suggests that our definition of satellites is robust. Encourag- 
ingly, the faint end of the satellite luminosity function in the 
high resolution simul ation matches that of the newly discov- 
ered SDSS satellites (|Koposov et al-lfeoOct ) remarkably well. 
This provides further support for the 'vw' winds in which 
w w oc a and r)^ oc cr~ 2 . 

We again estimate the minimum resolved circular ve- 
locity of the subhaloes in the high-resolution simulation as 
described in Section 13.2.11 The minimum resolved value is 
6.1 km s , while the minimum value of the subhaloes that 
host a satellite is 12.1 km s - . Hence, all the subhaloes host- 
ing satellites are well resolved in the simulation. The mini- 
mum V max of the satellites in these models is set by re ion- 
ization |Okamoto et al.ll2008al ; lOkamoto fc Frenkll2009l ). 

In the lower panels of Fig. 1101 we compare the chemical 



Aq-A-4 has similar resolution to Aq-D-HR and Aq-A-1 has 
~ 100 times better resolution than Aq-D-HR. 
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Figure 10. Comparison of the low and high resolution simula- 
tions of the 'vw5<r model in the Aq-D halo. In all panels, the 
high resolution (Aq-D-HR) case is shown by solid lines and the 
low resolution (Aq-D) case by dotted lines. Upper-left: star for- 
mation histories of the central galaxies, as in Fig. [4] Upper-right: 
the satellite luminosity functions within 280 kpc, as in Fig. [5] The 
thin black dotted line shows the luminosity function of the low- 
resolution satellites that are matched to the bright (My ^ —11.2) 
high-resolution satellites. Lower-left: the iron abundance as a 
function of satellite luminosity, as in Fig. [6] Lower-right: oxygen- 
to-iron abundance as a function of the satellite luminosity, as in 
Fig. E] 



properties of the satellites in the two simulations. Both the 
iron abundance and the oxygen-to-iron ratios are consistent 
at the two resolutions. The high resolution simulation pre- 
dicts that the luminosity-metallicity relation for satellites 
should extend down to at least M v ~ — 7 (lower- left panel). 
There is a hint in the lower-right panel that the oxygen- 
to-iron abundance ratio in the simulated satellites increases 
slightly with metallicity. This trend is possibly inconsistent 
with the Local Group dwarf spheroidal data, in which the 
alpha-to-iron abundance ratios appear to be almost con- 
stant or slightly in creasing with decreasing metallicity (e.g. 
iTolstov et aUl2009h . 

Comparing matched pairs individually, we find that the 
stellar mass in the high-resolution satellites is, on average, 
30% higher than in their low-resolution counterparts. (This 
trend is opposite to that for the central galaxy where the 
high- resolution simulation forms a less massive galaxy.) The 
metallicity in the high-resolution satellites is, on average, 
only 8% higher than in the low-resolution simulation. 

Finally, in Fig. 1111 we show the star formation histo- 
ries of a selection of satellites in the high-resolution simu- 
lation. The top and middle panels show the star formation 
histories of the matched counterparts shown in panels (a) 
of Fig. [8] These histories are qualitatively similar to those 
seen in the low resolution simulation (Fig. [8}, indicating that 
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Figure 11. As (a) of Fig. [8] but for the 'vw5o-' model in the Aq- 
D-HR halo. The top and middle panels show star formation histo- 
ries of matched high-resolution counterparts of the low-resolution 
satellites shown in (a) of Fig. \E\ The bottom panels show star 
formation histories of satellites that are not resolved in Aq-D. 



the time variability is not a resolution effect but a result of 
the winds. The sporadic nature of the star formation makes 
the effective star formation timescale longer (and lowers the 
alpha-to-iron ratio) in the 'vw' models compared to the 'cw' 
models. As we saw earlier, the star formation histories of 
the surviving satellites retain no obvious imprint of reion- 
ization. Even the satellites in the bottom panels of Fig. 1111 
which were not resolved in the low-resolution counterpart 
(Aq-D), continue forming stars past the reionization epoch. 
We note, however, that the gas content in subhaloes that has 
failed t o form stars was evapora ted at the epoch of reion- 
ization l|Okamoto fc Frenkll2009l ). This process is likely to 
be responsible for the suppression of cosmic star formation 
rate at the epoch of reionizat i on fo u nd by larger cosmo log- 
ical simulations (|Crain et al.1 120091 ; ISchave et al.1 l2010h - It 
is interesting to note that star formation in these smallest 
satellites (Kn ax < 20 km s _1 ) ceases at the time of accre- 
tion or earlier. We discuss in detail these sorts of effects on 
surviving satellites and failed satellites ( i.e. subhaloes that 
do no t host stars) in a companion paper |Okamoto fc Frenkl 
I2009T ). 

In summary, our high resolution resimulation of the 
'vw5ct' model shows that the statistical properties of the 



satellite population in our low resolution calculations are 
robust to a change in mass resolution of a factor of ten. 
Furthermore, the high resolution simulation allows us to ex- 
tend the model luminosity function by 4 magnitudes, thus 
enabling a comparison with recent data. 



4 SUMMARY AND DISCUSSION 

We have performed cosmological hydro-chemodynamic sim- 
ulations of galaxy formation in 3 Milky Way-sized ha loes 
taken from the Aquarius project (|Springel et al.l [20081 ). A 
primary aim of our study has been to constrain the nature 
of some of the feedback processes that must have operated 
during galaxy formation. We have done this by focusing on 
satellite galaxies whose shallow potential wells make them 
particularly sensitive to feedback effects. As a byproduct, we 
have been able to explore some of the controversial proper- 
ties of satellite galaxies, such as their luminosity function, in 
the context of the ACDM cosmology which we have assumed 
in our simulations. 

Although our simulations are amongst the largest of 
their kind performed to date, they lack the resolution to 
follow the physics of the interstellar medium directly. Such 
processes must be included as 'sub-grid physics'. In particu- 
lar, we have assumed that energy injected by SNe generate 
galactic winds. We have assumed further that all the energy 
released by SNe is deposited as kinetic energy in the winds. 

Early studies, as well as more recent ones have shown 
that many observed galaxy properties cannot be reproduced 
unless most of the SN energy is used to blow gas out of galax- 
ies (e.g. ISpringel fc H ern quistj l2003bl; lOkamoto et al.ll2005l; 
Oppenheimer fc Dave|[2008l ; ICrain et al. 2009; S chave et al.l 
20101 ) . Since a large fraction of SN energy must be radi- 



ated away, other sources of energy seem to be required 
to power strong winds such as radiatio n from young stars 
l|Murrav et all 120051 : iNath fc Silkl I2009T ). In this study we 
have included these unresolved processes in a purely phe- 
nomenological way by assuming two types of 'energy- 
conserving' winds. In the first type, the initial wind speed 
is taken to be proportional to the local velocity dispersion 
of the dark m atter, u w cx a, as suggested by recent data 
l|Martinll2005l ), and thus, the wind mass-loading, rj w cc a 
('vw' models). In the second type, the initial wind speed and 
the wind mass-loading are assumed to be constant for all 
galax ies ('cw' models), as suggested by earlier data l|Martinl 
Il999l ) . Disc-dominated galaxies formed in several of our sim- 
ulations. 

In some of our simulations, we included a multiphase 
model for star-forming gas, but we foun d that this made lit - 
tle difference to the outcome (see also ISchave et al.l |2010| ) , 
apart from the overall morphology of the central galaxy. 
The use of a stiff equation of state for the ISM tends to 
promote the formation of disc-dominated galaxies by stabil- 
ising gaseous discs against gravitational instability. The key 
to the diversity of behaviours that we find is not the treat- 
ment of the ISM, but rather the treatment of the SN-driven 
winds. 

The 'vw' models give a reasonable match to the ob- 
served luminosity function of Local Group satellites. A ma- 
jor factor in this success is the behaviour of the mass-loading 
in the wind which becomes increasingly large in smaller 
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galaxies. If the mass-loading is kept constant as in the 'cw' 
case, galaxy formation is not sufficiently suppressed in small 
haloes resulting in a satellite luminosity function which rises 
much too steeply at the faint end. Previous SPH simula- 
tions had already shown that an acceptable match to the 
abundance of brigh t satellites can be obta i ned from ACDM 
initia l conditions (iLibeskind et alj 120071 ; iGovernato et al.1 
120071 : iMaccio' et aUbOOSh . thus confirming the conclusion 
from early semi-analytic studies that the so-called ' satel- 
lite problem' in t he CDM cosmology (|Klypin et al.l 1 19991 ; 
iMoore et all 19991 ) disappears when proper account is taken 
of the baryonic processes involved in galaxy formatio n 
dBullock et al . 2000:; iBenson et al.ll2002j ; ISomervilld [2002). 
Our high-resolution simulation, however, extends the com- 
parison with the luminosity function to much fainter mag- 
nitudes than previous SPH simulations. 

The 'vw' feedback model that produces an acceptable 
satellite luminosity function, also produces an acceptable 
luminosity-metallicity relation for the satellites. By contrast, 
the 'cw' model fails to reproduce this relation for the same 
reasons that it fails to reproduce the luminosity function: 
the constant mass-loading leads to an iron abundance in 
satellites with v csc S> v w which is essentially independent of 
luminosity, contrary to what is observed. 

None of our 'vw' models reproduces the oxygen-to-iron 
abundance ratio measured for dwarf satellites. There are 
several possible reasons for this discrepancy: the star for- 
mation timescale may be too long, the model of SNe la or 
the IMF that we have adopted may be incorrect, or our 
assumed oxygen yield may be too low. Given the large un- 
certainties in the value of the yield, the IMF and SNe la 
rates, we do not consider this shortcoming of our models 
to be fatal although further work is required. For example, 
a top-heavy IMF in these metal-poor objects could provide 
the solution although it is unclear whether such a radical 
assumption would destroy the agreement of the model with 
other observables. 

Based on our estimate of the minimum resolved circu- 
lar velocity of subhaloes in the simulations, we find that 
subhaloes hosting satellites in the 'vw' models are reliably 
resolved down to the faintest magnitude. Thus, the faint-end 
slope of the satellite luminosity function and the luminosity- 
metallicity relation in the 'vw' models are robust. This con- 
clusion is confirmed by a direct convergence study of one 
model ('vw5<t' of the Aq-D): low- and high-resolution sim- 
ulations give consistent results. On the other hand, in the 
'cw' models in which the smaller wind mass-loading factor 
allows satellites to form in halos with small values of Vma.x, 
the faint-end slope of the satellite luminosity function could 
be affected by limited numerical resolution although we esti- 
mate that any such effects would be restricted to the extreme 
faint end. 

Barring the alpha-to-iron ratio, we conclude that an 
energy-conserving wind model in which the mass load- 
ing scales inversely with square of velocity dispersion pro- 
vides a viable model for feedback as judged by the proper- 
ties of satellite galaxies. Our preferred wind model is dif- 
fer ent from the 'momentum -driven' wind models favoured 
by iFinlator fc Dave] (|2008l ) in order to explain the mass- 
metallicity relation of large galaxies. In this kind of model, 
the wind speed is proportional to a as in our 'vw' mod- 
els, but the mass-loading scales as i|w oc it" 1 (rather than 



as 77 w <x a~ ). It could be that the mechanisms that drive 
galactic outflows in small and large galaxies are different. 

In 'vw' models with wind speed, v v ~ 4 — 5<r, the mass 
loading is large in small galaxies. Winds therefore remove 
substantial amounts of star-forming gas, but much of this 
material eventually falls back onto the galaxy. As a result, 
the star formation in small satellites is episodic and has a 
much longer timescale than in the 'cw' models in which the 
gas is expelled from the small halo, never to return. These 
general properties underlie the different predictions for the 
number density of faint galaxies and their chemical proper- 
ties in the different models. 

The star formation histories of the satellites retain no 
obvious imprint of the reionization of gas at early times. 
Accretion onto larger haloes, on the other hand, often af- 
fects subsequent star formation and, in low mass satellites 
(Knax < 20 km s _1 ), it truncates it altogether. Only 5% 
of subhaloes whose circular velocities are higher than the 
minimum resolved value (V ma x > 6.1 km s _1 ) in the high- 
resolution simulation (Aq-D-HR) host satellites. The vast 
majority of subhaloes do not manage to make a visible 
galaxy and remain dark. The mechanisms that distinguish 
visible from dark satellites are investigated in a companion 
paper (lOkamoto fc Frenkl [20091 ). 
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